#ifndef _HF_SOLVER_M_H_
#define _HF_SOLVER_M_H_

#include<iomanip>
#include<algorithm>//for sort algorithm
#include<cmath>
#include<eigen3/Eigen/Eigen>//matrix library

#include"orbital.h"
#include"system_m.h"
#include"help.h"

using Eigen::Matrix;
using Eigen::Dynamic;
using Eigen::SelfAdjointEigenSolver;

using Eigen::ComplexEigenSolver;
using std::cout;
using std::stable_sort;
using std::abs;
using std::setw;


///Hartree-Fock orbital, M 
class HF_Orbital_TzM:public Orbital_parmmtz
{
 public:
  typedef double DataType;
  HF_Orbital_TzM(){}
  HF_Orbital_TzM(int _par,int _mm,int _tz,DataType _e):Orbital_parmmtz(_par,_mm,_tz),e(_e){}

  void set(int _GroupIndex,int _par,int _mm,int _tz,DataType _e)
  {
    GroupIndex=_GroupIndex;
    par=_par;
    mm=_mm;
    tz=_tz;
    e=_e;
  }
  
  ///wave function under certain bases(eg. under harmonic bases).
  ///GroupIndex: index of the group in the groups of bases 
  ///coeff: coefficients under bases in that group
  int GroupIndex;
  vector<DataType> Coeff;
  bool state;///< 0 for hole, 1 for particle

  bool operator < (const HF_Orbital_TzM & Orb) const
  {
    //    return e<Orb.e;
    //in case that energy is too close
    if( abs(e-Orb.e) >1e-6)// numeric_limits<double>::epsilon() )
      {
    	return (e < Orb.e);
      }
    else
      {
    	return mm<Orb.mm;
      }
  }
  DataType e;///< single particle energy.
};




class HFSolver_M
{
 public:
  typedef System_TzM_HO System_MType;
  typedef typename System_MType::DataType DataType;
  typedef typename System_MType::OrbitalType OrbitalType;
  typedef Matrix<DataType,Dynamic,Dynamic> Mat;
  typedef Matrix<DataType,Dynamic,1> Vec;
  
  HFSolver_M(){}
  HFSolver_M(System_MType * Ptr):pSystem(Ptr){}
  HFSolver_M(System_MType * Ptr,int _A,int _Z):pSystem(Ptr),A(_A),Z(_Z){}
  void setup(System_MType * Ptr,int _A,int _Z)
  {
    pSystem=Ptr;
    A=_A;
    Z=_Z;
  }
  
  void setAZ(int _A,int _Z){A=_A;Z=_Z;}

  DataType getOneBody(int bra,int ket) const
  {
    return pSystem->get1B(A,bra,ket);
  }
  
  DataType getTwoBody(int a,int b,int c,int d) const
  {
    return pSystem->get2BV(A,a,b,c,d);
  }
  
  DataType HFDiagram(int bra,int ket) const;

  ///return interaction mean field after HF.
  DataType getVHF(int bra,int ket) const;

  
  ///single particle hamiltonian, <bra|V1B|ket> + 1/(2j+1) sum_h{ (2J+1) < bra h->J | V2B | ket h->J >}, h is hole
  ///(group,a) corresponding bra,(group,b) corresponding ket
  DataType HFhamil(int group,int a,int b) const
  {
    int bra=pSystem->Groups[group][a];
    int ket=pSystem->Groups[group][b];
    return getOneBody(bra,ket)+HFDiagram(bra,ket);
  }

  
  void printOrbitals(int num=-1) const
  {
    if(num==-1)
      num=HFOrbitals.size();
    cout<<"num"<<"\t"<<"par"<<"\t"<<"mm"<<"\t"<<"tz"<<"\t"<<"e"<<"\t"<<"OccuP"<<endl;
    for(int i=0;i<num;i++)
      {
      	cout<<i<<"\t"<<HFOrbitals[i].par<<"\t"<<HFOrbitals[i].mm<<"\t"<<HFOrbitals[i].tz<<"\t"<<HFOrbitals[i].e<<"\t"<<1-HFOrbitals[i].state<<endl;
      }
  }
  
  ///initialize
  void init();

  ///sort and fill orbitals.
  void fillOrbitals();

  ///construct matrix and diagonalize, generate new HF orbitals
  void diag();
  
  ///iteration
  void iter();

  ///cal. HF energy
  void calEnergy();
  
  int A,Z;
  vector< HF_Orbital_TzM > HFOrbitals;
  DataType Energy;

  System_MType * pSystem;
};

//*********************************************************************
#endif
